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THREE-DIMENSIONAL  NUMERICAL  SIMULATIONS 
OF  FELS  BY  TRANSVERSE  MODE  SPECTRAL  METHOD 

I.  INTRODUCTION 

Tliree-dimensional  simularions  of  the  FEL  gain  process  will  become  increasingly  im¬ 
portant  in  the  design,  optimization  and  interpretation  of  FEL  experiments.  Since  the 
fundamental  FEL  concepts  have  been  successfully  demonstrated,  the  FEL  is  noAv  pushing 
the  limits  of  the  FEL  capabilities  in  all  directions.  Due  to  the  increasing  complexity  of 
FEL  experiments,  simple  one-dimensional  theory  is  often  not  adequate.  At  the  same  time, 
however,  the  analytical  three-dimensional  results  are  only  obtainable  for  a  limited  number 
of  special  operating  conditions.  The  importance  of  the  three-dimensional  effects  will  vary 
with  the  experiments.  Some  examples  of  the  three-dimensional  effects  that  need  to  be 
understood  are: 

(a)  The  transverse  electron  beam  profile  is  often  asymmetric,  resulting  in  the  asymmetry 
of  the  radiation  field.  This  is  most  pronounced  in  the  storage  rings  and  circular 
microtrons, 

(b)  The  radiation  field  varies  both  in  the  transverse  and  axial  directions.  There  is  a 
crucial  trade-off  of  large  filling  factors  with  short  Rayleigh  lengths  and  vise  versa. 
The  optimal  situation  varies  with  the  experiments. 

(c)  Betatron  oscillations  from  the  wiggler  field  cause  electrons  to  sample  a  varying  radia¬ 
tion  field.  Under  certain  situations,  this  can  cause  betatron-synchrotron  instabilities. 

(d)  FELs  are  pushing  toward  high  gain  operation.  Strong  self-focusing  in  a  high  gain 
FEL  can  substantially  alter  the  wave  front  properties.  A  resonator,  not  self-consistly 
designed  with  FEL  physics,  most  liltely  is  not  the  optimal  design  for  the  high  gain 
FEL  operation. 

(e)  Sideband  instabilities  in  an  oscillator  can  lead  to  pulse  breakup  and  affect  the  qualiu* 
of  the  radiation  beam. 

The  analytical  and  numerical  methods  employed  in  the  study  of  the  transverse  \’aria- 
tion  of  the  wave  equation  have  taken  the  approaches  of  transverse  mode  spectral  method. 
Manuscript  approved  September  18, 1983. 


transform  spectral  method,®'^  spectral  method,"  finite  difference  method,®’^  and  Lienard- 
Wiechert  potential  method^*’.  The  different  methods  refer  to  the  techniques  used  to  evalu¬ 
ate  the  A/j  term  in  the  wave  equation,  where  A/i  is  the  vector  potential  of  the  radiation 
field.  All  three  kinds  of  spectral  methods  involve  representing  the  solution  of  the  radiation 
field  as  a  truncated  series  of  known  functions  of  the  independent  variables. 

The  transverse  mode  spectral  method  decomposes  the  radiation  field  into  a  truncated 
series  of  a  complete  set  of  orthogonal  functions,  which  satisfy  the  free  space  wave  equations 
with  the  appropriate  boundarj'  conditions.  One  takes  advantage  of  the  transverse  mode 
properties  to  reduce  the  wave  equation  to  a  set  of  first  order  differential  equations. 

For  the  spectral  method  and  transverse  spectral  method,  the  series  takes  on  the  form  of 
a  transform,  and  the  wave  equation  reduces  to  a  simple  form  in  the  transformed  variables, 
i.e.,  has  an  analytical  form  in  the  transform  space.  The  transform  spectral  method 

requires  numerically  taking  the  transform  of  the  driving  current.  The  spectral  methods 
describe  the  current  in  terms  of  the  Lagrangian  variables,  and  the  current  term  can  be 
evaluated  analytically. 

In  Sec.  n.  we  will  outline  the  transverse  mode  spectral  method  in  general.  i.e..  not 
specifying  the  form  of  the  transverse  modes.  This  formulation  will  include  slow  transverse 
motion  and  betatron  oscillations  of  the  electrons  in  the  realistic  wiggler,  finite  emittance, 
energy  spread  and  self-consistent  a.xial  particle  dynamics.  This  method  conserves  energy. 
The  formalism  is  extended  to  the  study  of  sideband  formation^"^’^^"^"  on  a  long  electron 
pulse  in  the  FEL  oscillator. 

In  Sec.  HI,  we  outline  the  spectral  and  the  transverse  spectral  methods.  The  three 
different  spectral  methods  are  compared  in  this  section. 

The  ad'vantages  of  the  transverse  mode  spectral  method  are;  i)  Free  space  wave  prop¬ 
agation.  finite  size  mirrors,  and  apertures  can  be  handled  analytically,  ii)  It  is  easy  to 
include  transverse  particle  motion  exactly,  iii)  This  method  lends  itself  to  analytical  and 
semi-analytical  solutions  and  can  provide  physical  insight  for  many  problems,  iv)  The 
transverse  boundary  conditions  are  included  automatically  in  the  waveguide  mode  expan¬ 
sion. 
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In  Sec.  IV,  we  apply  the  transverse  mode  spectral  method  to  an  example  with 
self-focusing  properties.  In  this  example  the  electron  beam  focuses  due  to  gain  and 
refraction. 

n.  THE  TRANSVERSE  MODE  SPECTR,\L  METHOD 
A.  THE  WAVE  EQUATION 

In  this  paper,  we  will  consider  only  the  linearly  polarized  wiggler,  since  that  is  the 
most  common  wiggler  field.  The  formulation  for  a  circularly  polarized  wiggler  requires 
only  minor  modifications.  The  realistic  magnetic  field  will  be  expressed  in  terms  of  the 
vector  potential  of  the  wiggler, 


,(y, ::)  =  Aa,iy,=)cos{  /  *,„(5')d-')ex, 

Jo 


where  :r)  =  .4ii/(r) cosh(/:u^y),  and  are  the  slowly  varying  amplitude 

and  wavenumber  of  the  wiggler. 

W’hen  kyjij  <<  1.  one  can  use  the  approximation  .4^,  -4u((-)  for  the  calculation 

of  gain  and  axial  particle  dynamics.  The  accurate  form  for  .4u/  is  only  necessary  in  the 
calculation  of  the  betatron  oscillations.  For  convenience,  we  define  a  dimensionless  wiggler 
parameter  IT ( r)  =  {\e\  '  ^moc'^)A^(:). 

We  also  include  a  DC  accelerating  electric  field. for  the  purpose  of  efficiency 
enchancement.  For  convenience,  we  define  a  dimensionless  parameter  e£)c  =  l^l-E^uc/z^oC". 

The  radiation  field  will  be  expressed  as 


Aff{z,  y.  z,t)  =  -A(z,y,  -) 


exp[*(&j  - 


where 


L  \{ 


Aiz.y. ,')  =  ,4j(z.y,  ;)e.  ^  A,j{£,,!j,  z)ky  =  T^Y^  A.:...„(-)ee,.v,(z,  y, 

isO  m=:0 

~  "i" 

is  the  general  form  of  the  transverse  mode.  .4^,^  =  |-4t,mi«"’'-'"  is  the  complex  amplitude  of 
the  normal  modes,  Gi^m  =  represents  the  x-component  of  the  complex  transverse 

mode.  F/>,m  is  the  y-component  of  the  comple.x  transverse  mode,  (which  may  be  zero.) 
Aj(z.y.  ;)  =  ].4i(z.  y.  r)|  exp(j^(z.  y.  r)).  and  k  =  ^-/c.  We  use  lower  case  to  denote  the 
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normalized  parameters  of  the  radiation  field,  such  as  ax{x.  y,  z]  =  {\e\l \/2moC^)Ax(x,  y,  z) 
and  =  [\e\l yj2moc'^]At,m. 

The  function  should  satisfy’  the  free  space  wave  equation  and  the  appro¬ 

priate  boundary  conditions,  i.e., 

hi  +  ^  ^  0.  (3) 

The  orthogonality  condition  is 


J  j  ^t,m  '  ■< 


where  the  integration  is  over  the  appropriate  radiation  domain. 

The  wave  equation  for  the  FEL  is 

Substituting  the  expression  of  the  radiation  field  into  the  wave  equation,  and  making  the 
slowly  varying  approximation,  we  obtain 

V  +  C.C.  =  (6) 

,  oz  c 

t.m 

Multiply  both  sides  of  Eq.  (6)  by  exp(— ifAr  -  jjt))  and  integrate  in  time  over  a  period. 
Then,  dot  both  sides  by  and  integrate  over  the  transverse  dimension.  The  result  is 


^At_rr 

dz 


i  a-  /'  /* 


Now  we  need  to  e'valuate  the  current.  The  fiuid-like  electron  density  is 
D(x.!J.  :.Pz.Py.pZ-t}  =  <io{*0,yo-.Px.O-Py,0,Pz.O,  'r'ojHj!  - 

<^(Pl  -Pl)<5(Py  -P!,)(5(P,-  -P.-)-  (S) 

The  variables  with  and  the  variable  r  are  functions  of  the  Lagrangian  variables 
•fo- yo-Px.o- Pv.o.Pr.o,  -•  where  Xo^ijo  are  the  initial  transverse  positions  at  the  entrance 
of  the  wiggler  r  =  0,  Px.o^Py.o  are  the  initial  transverse  momentums  at  ;  =  0.  p^.o  is  the 
initial  a.xial  momentum  spread  at  r  =  0,  Wg  =  —jjto  is  the  initial  phase  of  the  electron  in 


the  ponderomotive  potential  well,  is  the  time  the  electrons  enter  the  wiggler  at  r  =  0, 
and  r  is  the  time  that  the  electron  reaches  position  z.  The  function  do  is  the  initial  phase 
space  density  distribution,  including  the  effects  of  emittance  and  energy  spread. 

The  Suid-like  beam  density  can  be  defined  as 


n{x.y,:J.}=  I  dp:.  I  dp,j  I  dp;D{x.  y.  :.p2.  p,j.  p,J). 

J  2C  *  ^  ^ 


We  will  define  the  effective  area  of  the  electron  beam  to  be 


<Tb{zJ)  =  ^  J  dx  j  dy  n{x,y,z,t), 

where  rio  =  n{x,y.  z,  t)\ma.x- 

The  current  density  can  be  defined  as 

=  -|e|i'.-o  J  dpi  J  dpy  J  dp:D{x.y,z,p2,Py.p:J)^^ 

=  -\e\noVzo  j  dto  I  dxo  j  dijo  J  dpx^o  J  dpy^o  J  dp,_o 

px  1 

- do{xo,yo-px.o^Py.o^Pz.o-  ii’o)d(x  -  xjd'(y  -  y)d(f  -  rj. 

pz  rio 

The  final  form  of  the  wave  equation  is 


—  =  ,(7F,(.-)A(.-|^ - ^ - 


where 


dwo  r ,  r  X  r  ^  r  ^  r  ^ 

<  ■">=  —  /  ^P--°  / 

Jo  ““  J —<x  J —CO  •'—00  J —CO  Jo 

(...)  do[Xo.  yo' Pz,0‘ Py,o- Pz.oi  '^0  )  (1«>) 

no<Tb 

is  the  ensemble  average,  the  normalization  is  <  (1)  >=  1.  C  =  ^’6  = 

(driei'^no/mo)^/'.  Fi(r)  =  Jo(b]  -  Ji(b)  comes  from  the  fast  electron  oscillation  in  the 
z-direction  due  to  the  linearly  polarized  wigjler,  6(;)  =  A''/2(l  -f  A'*),  w  =  JJ(A\.(z')  + 
A  -  -r  t'o,  d-  =  b~/c'.  Jzo  =  ’'zo/<^-  and  I’^o  is  the  mean  initial  a.\ial  velocity. 


B.  P.4IITICLE  DYN.AJvnCS 

The  mean  transverse  particle  position  and  momentum,  in  general,  should  be  integrated 
along  with  the  axial  particle  momentum  equations.  Under  certain  circumstances,  analytical 
forms  for  the  transverse  particle  trajectories  are  sufficient. 

Talcing  the  betatron  oscillations  as  an  etcample,  the  particle  motion  in  the  .x-direction. 
without  any  additional  external  focusing,  can  be  written  as 

x  =  Xo  +  {pxollomoc)z.  (14) 

The  particle  motion  in  the  y-direction  exhibits  betatron  oscillations.  Taking  «  1. 
the  particle  dj-namics  in  the  y-direction^®^  is  approximately 

where  kj  =  ijs  =  {pl,o^UzV^~- Py.o  =  Py,o ■  and  03  =  cos'M.'/c/.^/f)- 

The  particle  motion  in  the  z-direction  is  best  written  in  terms  of  the  equation  for  the 
total  relativistic  gamma  factor 

^  =  -TFi(r)/v(r)|a-(5..y.  r)|sin(j/.-^^(i.i/.  r))  -  (IS) 

dz  7 

and  the  equation  governing  the  phase 

(fw  (1  +  3zo)K  ^  ^ 

JU  = - i - (— (10 


where 


Wdz  ^  T 


is  the  degree  of  taper  for  the  efficiency  enhancement  schemes. 

47^  ki(z\  dz 

includes  the  etfect  resulting  from  the  combination  of  betatron  oscillation  and  contoured 


field,  and 


-5(.7o.Pyo.--)  =  2  /  fcd(-')<i--'+2o,. 

Jo 


'X'} 


6 


Equations  (12,  16-17)  form  the  self-consistent  set  of  equations  for  the  radiation  field  and 
the  particle  dynamics. 


C.  CONSERVATION  OF  ENTRGY 


We  can  show  that  this  formulation  conserves  energy'.  To  do  this,  we  rewrite  Eq.  (16) 
in  terms  of  the  amplitude  of  each  mode. 


dz 


-,F,(z)K{=\ 


X) sin  -eoc. 

t,m 


(20) 


where  is  the  phase  associated  with  each  mode.  VTe  have  assumed 

here  that  Re{,3i,m)^  <<  2t.  The  spatial  rate  of  change  of  the  total  electron  beam  energy 
is 


Lm 

~~  {moc-eoc) , 


(21) 


where  D  —  (<7ft/St:)(u,f /c'),i;oA.  The  spatial  rate  of  change  in  the  radiation  energy  is 


d:  Sjt  ^  dz  ^ 

t^m 


Applying  Eq.  (12)  to  Eq.  (22),  we  can  show  that  the  energy  is  conserved,  i.e.,  the  energy 
gained  by  the  radiation  is  equal  to  the  energy  lost  by  the  electrons  plus  the  energy  supplied 
by  the  axial  DC  electric  field. 


dlR 

dz 


(23) 


D.  THE  CHOICE  OF  THE  TR.AN'SVERSE  MODES  OF  THE  WAVE  EQUATION 

The  choice  of  the  normal  modes,  however,  will  depend  on  the  geometiy  of  the  FEL  of 
interest.  A  few  examples  of  the  possible  choices  of  the  normal  eigenmodes  are  given  below; 
(a)  Gaussian-Hermite  expansion  allows  the  most  flexibility  in  the  modeling  of  the  FEL 
physics  of  interest  in  an  open  resonator.  The  e.xpressions  for  are 
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Qt,m  — 


u  j)> 


(24a) 


't.m 


=  +r)s‘  -ii-tm  ^Ijtan 


-1 


(246) 


where  x  =  v'2x/'t’(-!.  j  =  \/2y/<t’(r),  Ht  is  ihe  Hermite  poiynomiai  of  :he  /.th  order. 


;  =  -  Zc)I^Q-  -c  is  tlie  axial  location  of  the  minimum  waist.  a'{r)  =  a'o(l  -i-  ■ 

Wo  is  the  spot  size,  and  Zq  =  w~k(2  \a  the  Rayleigh  length. 

(b)  When  the  electron  beam  and  the  resonator  system  are  axially  symmetrical,  and  when 
the  betatron  oscillations  are  not  included  in  the  model  of  the  linearly  polarized  wiggler, 
then  the  radiation  can  be  assumed  to  be  axially  symmetrical.  The  Gaussian-Leguerre 
modes  are  the  choice  for  the  expansion,  and  each  term  is  expressed  as  eo,m  =  Go,m^x, 
where 


Qo.m 


)exp(- 


(25a) 


^o.m  = -(2m -i- 1)  tan"^(i')  + (256) 
is  the  Laguerre  polynomial. 

(c)  In  a  low  gain  oscillator,  the  best  choice  would  be  the  eigenmodes  of  the  resonator 
cavit>'^®  including  the  effects  of  apertures,  finite  size  mirrors  and  resonator  losses.  The 
eigen  modes  in  turn  can  be  written  in  terms  of  the  appropriate  Gaussian-Laguerre  or 
Gaussian-Hermite  modes.  The  advantage  of  going  to  the  resonator  c.avity  modes  is 
that  the  losses  for  the  higher  order  resonator  cavity  modes  increase  rapidly  with  the 
cavity  mode  number,  such  that  only  a  few  modes  need  to  be  kept  in  the  caluclation. 
In  addition  there  is  no  need  to  calculate  the  radiation  outside  the  FEL  wiggler. 

(d)  Closed  waveguides  can  be  used  to  confine  the  radiation  beam  over  a  distance  long 
compared  to  a  Rayleigh  length.  If  the  waveguide  is  not  rectangular,  the  vector  poten¬ 
tial  of  the  wave-guide  modes  generally  have  both  .x  and  y  components.  For  the  analysis 
presented  here,  we  have  to  assume  the  transverse  guide  dimensions  are  large  compared 
to  a  radiation  wavelength.  The  theory  is  also  restricted  to  low-order,  low-loss  modes 
whose  propagation  constants  k  -f-  i?e(.i*..7,)  are  nearly  equal  to  the  plane-wave  value, 
i.e.,  E.e[3i^rn)  «  k.  For  wavelengths  that  are  comparable  to  the  waveguide  diameters 
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or  when  multimode  analysis  are  desired,  this  method  has  to  be  modified  slightly  to 
include  the  different  frequencies  and  axial  wavenumbers  associated  with  the  different 
waveguide  modes.  The  formulation  should  be  similar  to  that  described  in  the  Sec. 
ILF  with  sidebands. 

E.  COMPUTATION  SPEED 

The  cost  of  computation  per  axial  step  is  spent  mainly  for  the  radiation  field.  The 
number  of  operations  is  0!i.Ye(i  -r  llf-V  +  1)  to  integrate  Eq.  (12),  where  is  the  total 
number  of  electrons,  and  Oi  is  a  numerical  constant.  The  particle  orbit  equations  require 
ao(I  +  1){-V  -r  l).Y,/:Yo  number  of  operations  to  evaluate  the  radiation  field  a^lx.y.  z]. 
and  a^Nf.  number  of  operations  for  the  integration,  where  No  is  the  number  of  electrons 
per  ponderomotive  wave,  and  ao  and  a*  are  numerical  constants.  In  general,  Eq.  (16)  is 
faster  than  Eq.  (20).  which  requires  0'4(Z  1)(-V  —  l).Ye  number  of  operations. 


F.  freqtjtncy  sideb.myds  and  pulse  slippage 

The  electron  oscillation  in  the  ponderomotive  potential  well  is  called  the  synchrotron 
oscillation.  We  define  A'jo  to  be  the  synchrotron  wavenumber  of  the  electron  traveling 
exactlv  along  the  z-axis 


7-  f  ^k,yk  IV o  I 

A,o  =  - - —  A(r  a, 

\  7o 


When  the  intensiu'  of  the  radiation  becomes  high,  such  that  the  electrons  make  half  a 
bounce  in  the  wiggier.  then  sideband  frequencies  can  grow.  One-dimensional  simulations^" 
and  three-dimensional  simulations^"'  show  that  the  amplitude  of  the  radiation  field  be¬ 
comes  chaotic,  and  the  quality  of  the  radiation  field  becomes  degraded. 

In  order  to  understand  the  effect  of  sidebands  on  wavefront  curvatures,  the  three- 
dimensional  formulation  is  necessary.  We  assume  periodic  boundary,’  conditions  for  the 
radiation  field  of  length  U;,.  which  is  chosen  to  be  many  times  the  ampiitude  modulation 
distance  and  the  pulse  slippage  dbtance.  If  the  length  of  the  electron  beam  is  much  longer 


than  £e6,  only  a  section  of  a  long  electron  beam  needs  to  be  modeled.  If  the  electron  beam 
is  shorter  than  itb,  then  the  whole  pulse  shape  must  be  included. 


The  radiation  field  can  be  written  as 


Aaij.//.  .-./)  =  -A{x.!j,  ^J) 


exp  [iki] 


where 


A(x,»/.  ^.t)  =  z)e,  +  r)ey 


V  L  .V/ 


n=0  /=<)  /7i=0 

is  a  slowly  varying  function  of  position  following  the  radiation  pulse,  k  =  (1  +  i-o  )(7o /  (I  + 
K‘)]k^  is  the  resonant  wavenumber,  Akn  =  {n  -  l)6k,  5k  =  and  ^  =  r  -  ct. 

A(x.  !j.  f)  is  3  slowly  varj-ing  function  of  position  and  time  following  the  radiation  pulse. 

Following  the  same  procedure  outlined  in  Sec.  II.A  and  II.B.  we  obtain  the  self- 
consistent  set  of  equations  for  7  and  w.  The  independent  Lagrangian  variables  are 

^0'  !Jo  •  P2.0  •  Py.Q  •  Pz,o^  ^0 1  ^0"  Und 


c  dt 


^iCF,{:]K[z] 


-37  =  -T-ri(-)A*!aj(i..y.  0| sinfji' +  y,  0)  - 

c  dt  7 


where 


dt^  7  c  dt 


=  f  ^  ^  I  I  I  I  ^P'-° 

Jo  Jo  J  J  J — ro  J  —  yz  Jo 


(...)  d(j{Xo,  ^Q.  Pj^q,  Py^Q.  P;^J,  li'o’  ^0]  (31) 

tlo^b 

is  the  ensemble  average.  =  [2'xj(lib]so-:  ^Jid  0  <  <  l.b  are  the  axial  location  of  the 

electron  in  the  electron  pulse  at  t=0,  ip  =  jl[k^[ct')+k]vzdt'  -  ckt  +  iPo,  z  =  /J  +  ^o. 


^  —  ct  ^ol5k  —  s,  and  s  =  ef(l  +  K‘)f2'';l  is  the  pulse  slippage  distance. 

The  method  described  above  is  straightforward.  Since  the  number  of  electrons  that 
must  be  used  to  model  the  electron  beam  with  finite  emittance  and  finite  length  is  large, 
any  general  numerical  method,  including  this  method,  requires  a  considerable  amount  of 
computation  time.  For  limited  computational  resources,  it  is  necessary  to  form  simplified 
models.  This  formulation  lends  itself  towards  semi-analytical  models  of  finite  length  elec¬ 
tron  pulses/“^’^®  where  some  of  the  integrals  in  the  expression  for  ensemble  average  can 
be  evaluated  analytically. 

m.  THE  SPECTRAL  .\ND  TR.4NSFORM  SPECTRAL  METHODS 

The  principle  behind  the  spectral  and  transform  spectral  methods  is  similar  to  the 
transverse  mode  spectral  method.  The  major  difference  comes  from  the  representation 
of  the  radiation  field,  for  example.  Fourier  series  in  the  Cartesian  geometry  and  Hankel 
transform  in  the  axially  cylindrical  geometry. 

Let  us  take  the  Fourier  series  for  the  illustration.  The  domain  for  the  radiation  field 
is  -Dxl2  <  X  <  Dxi2  and  -Dyj2  <  y  <  Dy/2,  and  the  boundaiy-  conditions  are  periodic. 
The  slowly  varying  amplitude  of  the  radiation  field  is  written  as 

.  L  M 

A(a:,j,.-)  =  E  (32) 

«=0  m=0 

where  ki  =  2zijDx  and  km  =  2:rmjDy.  The  equation  for  the  amplitude  in  terms  of  the 
current  is 


When  the  current  Jx  is  formulated  in  Eulerian  variables  for  the  x  and  y  coordiantes. 
then  the  right-hand-side  of  Eq.  (33)  has  to  be  evaluated  numerically,  and  th<^  method  is 
called  the  transform  spectral  method.  If  the  current  is  formulated  in  Lagrangian  variables 
(^o-yo-Px,o-Py,o-Pz.j’  then  Eq.  (33)  can  be  reduced  to 
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d  A?  it'  1  J  p-i[k,i+ii„y-<u)\ 

T.  +  '^J  =  ••<7n(--)A-(.-)  - ) 


Tlie  particle  dynamics  equations  (16-17)  are  also  applicable  here.  For  a  Fourier  series 
expansion,  one  can  show  that  the  enerjj*  is  conser\'ed  in  the  D.  by  D.j  domain  for  the 
spectral  method. 

In  general,  the  Eulerian  formulation  of  the  current  is  inferior  to  the  Lagrangian  for¬ 
mulation.  The  Eulerian  formulation  is  not  convenient  to  use  for  the  study  of  betatron 
effects,  because  the  numerical  transforms  of  Eq.  (33)  require  the  current  to  be  evaluated 
at  prespecified  grids.  This  problem  can  be  reduced,  if  the  number  of  grids  used  across  the 
electron  beam  is  large,  which  in  turn  requires  a  larger  number  of  terms  in  the  expansion. 
For  the  Hankel  transform,  there  is  the  additional  problem  of  numerical  errors  resulting 
from  numerical  integration  with  a  finite  number  of  grids  in  an  infinite  integration  domain. 

The  estimate  for  the  speed  of  computation  for  the  spectral  method  has  the  same  form 
as  in  Sec.  II.E.  The  value  for  the  coefficient  an  and  a*  depends  on  the  expansion,  and 
they  are  smaller  for  the  Fourier  transform  and  larger  for  the  Hankel  transform. 

The  transverse  mode  spectral  method,  in  some  instances,  is  better  than  both  the  spec¬ 
tral  and  transform  spectral  methods,  because  of  the  following  properties.  The  propagation 
of  the  radiation  field  through  apertures,  and  the  refiection  and  transmission  of  the  radia¬ 
tion  field  at  mirrors  can  be  evaluated  in  terms  of  matrix  multiplication  of  the  amplitudes  of 
the  normal  modes.  The  spectral  and  transform  spectral  methods  require  an  additional  step 
in  converting  the  radiation  field  into  Gaussian  modes.  Finally,  the  transverse  mode  spec¬ 
tral  method  is  easier  to  implement  for  calculations  in  complicated  waveguide  geometries, 
because  the  boundary  conditions  are  automatically  included. 

The  transverse  mode  spectral  method  is  not  the  best  numerical  scheme  when  the 
FEL  wiggler  is  many  times  the  Rayleigh  length  and  the  FEL  radiation  is  strongly  focused 
to  the  electron  beam.  The  reason  is  that  the  higher  order  modes  become  more  and  more 
important.  The  number  of  transverse  modes,  which  have  to  be  included,  become  large.  The 
Fourier  series  spectral  method  employing  the  Lagrangian  Formulation  is  more  appropriate 
for  this  type  of  simulations. 


•  •  •  •  -  •  -O  ^ 
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IV.  SELF-FOCUSING  EX.A_\IPLE  .-LN'D  DISCUSSION 

A  three-dimensional  code  employing  Gaussian-Hermire  expansion  of  the  radiation  field 
is  applied  to  a  linearly  polarized  higgler.  The  following  examples  will  show  self-focusing 
properties  of  the  FELs. 

The  electron  beam  is  assumed  to  have  a  Gaussian  profile  and  a  radius  of  2.25  x  10~^ 
cm.  The  current  is  50  A  and  the  energy  is  109.6  MeV.  The  wiggler  has  a  magnetic  field 
5a,  =  6.3  kG,  period  of  2.4  cm  and  length  !„,  =  6  m.  The  incident  0.5  nm  radiation  at 
the  entrance  of  the  wiggler  is  a  Gaussian  TE\Iqq  mode  with  a  minimum  waist  Wq  =  0.06 
cm  located  at  r  =  Ljujl.  and  a  power  of  5.8  x  10^  W. 

In  this  example,  the  electron  beam  radius  is  much  smaller  than  the  spot  size,  i.e., 
Teb  «  Wo-  If  the  radiation  has  the  spot  size  of  the  electron  beam,  its  Rayleigh  length 
would  be  only  32  cm.  A  plot  of  the  power  gain  is  shown  in  Fig.  1.  The  FEL  is  operating 
in  the  high  gain  regime,  and  the  radiation  saturates  at  4  m. 

Figures  2.a-e  shows  the  radiation  amplitude  at  z=0.  150,  300,  450  and  600  cm  plotted 
from  -4wo  to  4u;o  in  both  the  x  and  y  plane.  In  this  example,  the  self-focusing  phenomenon 
is  not  only  due  to  refraction,  but  also  gain.  The  peak  amplitude  in  Figs.  2a-e  are  4.5x10“^. 
1.1  X 10“^,  4.6  X 10“^,  4.1  X 10“^  and  2.6  x  10~^  respectively.  At  150  cm,  the  laser  beam  is 
already  narrowed  down  significantly.  The  corresponding  phase  front  is  shown  in  Fig.  3.a. 
The  focusing  due  to  the  electron  beam  is  manifested  by  the  small  mound  at  the  center. 

The  beam  radius  remains  at  roughly  the  same  size  from  150  cm  to  the  end  of  the 
wiggler.  Due  to  electrons  oscillating  away  from  the  bottom  of  the  ponderomotive  potential 
well,  the  radiation  loses  energ3.'’  4  m.  In  Refe.  13-14,  it  was  shown  analytically  that 
FEL  radiation  is  not  only  governed  by  gain  and  diffraction,  but  also  by  refraction  when 
operating  in  the  trapped  particle  mode.  For  this  case,  the  resonant  phase  is  at  the  origin, 
i.e.,  sin  il'h  =  0.  Even  though  the  radiation  is  losing  energy  after  r  =  4  m,  the  radiation 
beam  around  the  electrons  still  focuses,  while  the  radiation  further  away  from  the  electrons 
defocuses,  as  shown  in  Fig.  2.e.  This  b  in  qualitative  agreement  with  the  theory  given  by 
Refs.  13-14.  Figure  3.b  shows  the  focusing  properties  of  the  wave  front.  It  is  important 
that  the  radiation  beam  does  not  defocas  when  the  radiation  is  losing  energy  due  to  the 
bouncing  of  the  electrons  in  the  ponderomotive  potential  well. 
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The  FEL  radiation  does  not  always  undergo  focusing,  for  example,  the  electron  beam 
could  dig  a  hole  in  the  laser  amplitude  profile  when  the  FEL  is  operating  with  loss  in  the 
low  gain  regime  and  large  frequency  mismatch. 

Xe.xt.  we  give  an  example  where  the  center  of  the  electron  beam  does  not  travel  down 
the  center  of  the  wiggler  due  to  irregularities  in  the  magnetic  wiggler  field.  The  center  of 
the  electron  beam  is  assumed  to  execute  a  slow  sinusoidal  motion.  Ucntr  =  sin(rc/Iu^ ). 
The  radiation  guiding  properties  of  the  e-beam  is  still  clearly  evident.  Fig.  4.  The  peak 
amplitude  is  2.4  x  10“^.  The  motion  of  the  center  of  the  electron  beam,  however,  causes 
significant  distortion  of  the  radiation  beam. 

In  summary,  this  paper  outlines  one  method  of  solving  the  three-dimensional  FEL 
radiation  field  self-consistently  with  the  electron  dynamics.  The  advantages  of  this  method 
are:  1)  The  boundaries  in  the  transverse  directions  are  included  automatically.  2)  It  is  easy 
to  incorporate  transverse  particle  motion.  3)  Free  space  prapagation,  finite  size  mirrors 
and  apertures  can  be  handled  analytically.  4)  This  method  lends  itself  to  analytical  and 
semi-analytical  solutions. 

The  disadvantages  of  this  method  is  the  increase  in  the  computation  time  if  a  large 
number  of  electrons  are  used  in  the  radial  direction,  and  if  a  large  number  of  modes  are 
required. 

The  numerical  examples  illustrate  the  self-focusing  property  of  the  FEL.  Under  appro¬ 
priate  conditions,  the  laser  radiation  maintains  a  roughly  constant  radius  radiation  beam. 
We  also  showed  that  the  irregularities  in  the  wiggler  field  can  cause  significant  distortion 
of  the  radiation  beam.  The  implications  of  this  for  the  FEL  design  are  important.  For 
electron  beams  with  good  emictance.  it  is  possible  to  focus  the  beam  to  a  small  area  and 
operate  the  FEL  in  the  high  gain  regime,  .and  not  have  to  worr>’  about  difiiraction.  In  this 
case  it  is  important  to  design  the  resonator  cavity  including  the  FEL  gain  model. 
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Fig.  2  (Corn’d)  —  Plots  of  the  amplitude  of  the  radiation  field  in  the  x  and  y  plane  at 
(a)  z  -  0,  (b)  2  -  150  cm,  (c)  z  -  300  cm,  (d)  z  -  450  cm  and  (e)  z  -  600  cm. 
The  tick  marks  are  separated  by  a  distance  w,. 


Fig.  2  (Cont’d)  —  Plots  of  the  amplitude  of  the  radiation  field  in  the  x  and  y  plane  at 
(a)  z  —  0,  lb)  z  —  150  cm,  (c)  z  —  300  cm,  (d)  z  -  450  cm  and  (e)  z  —  600  cm. 
The  tick  marks  are  separated  by  a  distance  Wg- 


Fig.  2  (Cont’d)  —  Plots  of  the  amplitude  of  the  radiation  field  in  the  x  and  y  plane  at 
(a)  z  *•  0,  (b)  z  —  150  cm,  (c)  z  *  300  cm,  (d)  z  =-  450  cm  and  (e)  z  *  600  cm. 
The  tick  marks  are  separated  by  a  distance  w„. 


Fig.  2  (Corn’d)  —  Plots  of  the  amplitude  of  the  radiation  field  in  the  x  and  y  plane  at 
(a)  z  *  0,  (b)  z  ”  150  cm,  (c)  z  *  300  cm.  (d)  z  *  450  cm  and  (e)  z  *  600  cm. 
The  tick  marks  are  separated  by  a  distance  w,. 


Fig.  3  —  Plots  of  the  phase  of  the  radiation  field  in  the  x  and  y  plane  'a'  z  =  150  cm 
and  <b)  z  =  cOO  cm.  The  tick  marks  are  separated  by  a  distance  The  discontinuity 


is  an  artifact. 
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